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Qli| Abstract 

We compare the phase space slicing and dipole subtraction methods in the computa- 
tion of the inclusive and differential next-to-leading order cross sections for heavy quark 
production in the simple process 7* — > QQ. For the phase space slicing method we study 
■ the effects of improvement terms that remove restrictions on the slicing parameter s m i n . 

For the dipole method our comparison is a first check on some of its counterterms in- 
volving massive quarks, derived recently. In our comparison we address issues such as 
numerical accuracy and efficiency. 



1 Introduction 



Fully differential QCD cross sections are important observables for studies at high-energy col- 
liders. By allowing detector-specific acceptance cuts on phase space variables they eliminate the 
need for extrapolation into unmeasured, and often also poorly calculable regions, and thereby 
improve theory-experiment comparisons. Reliable theoretical predictions for such differential 
cross sections require the inclusion of at least next-to-leading order (NLO) QCD corrections. 
NLO calculations combine virtual one-loop corrections with the real emission contributions 
from unresolved partons. These two parts are usually computed separately and each is infrared 
divergent, only their sum is infrared finite. NLO Monte-Carlo programs incorporate both pieces 
and allow the simultaneous computation of many differential cross sections for the particular 
reaction considered. 

However, these programs require that infrared singularities be eliminated before any nu- 
merical integration can be done. There are essentially two types of methods to effect this 
cancellation. The phase space slicing (PSS) method [[I], 0, |3], f|] is based on approximating the 
matrix elements and the phase space integration measure in boundary regions of phase space 
so integration may be carried out analytically. The subtraction method |5], |7| is based on 
adding and subtracting counter terms designed to approximate the real emission amplitudes in 
the phase space boundary regions on the one hand, and to be integrable with respect to the 
momentum of an unresolved parton. 

For massless partons both methods are well-developed and have been widely used. A quite 
general formulation of phase space slicing has been given in Ref. [§,0- It was extended to 
include massive quarks and identified hadrons in Ref. [|l(| . 



There exist two general formulations of the subtraction method. One is the residue approach 



TT| , the other the dipole formalism [ 12| . Both can handle massless partons and identified 



hadrons in the final and/or initial state. The extension of the dipole method to handle massive 



quarks, using dimensional regularization, has been given recently in Ref. [13 . An extension 



to photon radiation off massive fermions, using small masses for infrared regularization, was 



developed by Dittmaier in [14]. There are also hybrid methods [15] that combine elements of 
the slicing and subtraction methods such that both the resolved and unresolved contributions 
are numerically small and can be reliably integrated. 

With general formulations of the phase space slicing and dipole methods for massless and 
massive quarks now available, it is interesting to compare their efficiency and accuracy. In this 
paper we do this for (differential) "cross sections" for heavy quark production in the process 
7* — > QQ. This case is of course very simple but also generic for more complicated processes. 



In the case of the NLO cross section for UH production fLq , |17|] it was recently verified |L7 
that the slicing method and a somewhat differently phrased dipole method |18| agreed. 

Our results using the dipole subtraction method represent the first numerical implemen- 
tation of some of the subtraction terms computed in While the dipole method is exact, 
the PSS method requires the introduction of a theoretical resolution parameter s min , usually 
required to be quite small. We include improvement terms in the PSS method [jn|, and study 
their effect of removing restrictions on the size of the phase space slicing cutoff. This is espe- 
cially important for cross sections involving heavy quark production, and allows for a free choice 
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of slicing parameter without reference to the heavy quark mass, a prerequisite for considering 
the high-energy or zero- mass limit. 

The paper is structured as follows. In section 2 we compute the fully differential cross 
section for 7* — > QQ using the PSS method, including the improvement terms. In section 3 
we compute this cross section with the dipole subtraction method. In section 4 we present a 
numerical comparison of the two methods, followed by our conclusions. 

2 Phase Space Slicing 

We consider the process 

l*(q)^Q(pi) + Q(p2), (1) 

with p\ = p\ = m 2 . The NLO corrections involve virtual corrections to ([I]) and the gluon 
bremsstrahlung reaction 

i*(q)^Q(pi) + Q(P2) + g(P3), (2) 

with p 2 = 0. We define the invariants 



s 



2pi-pj, = {pi+Pjf . (3) 



'j 



The final state phase space for the 3 parton contribution is divided into "hard" and "soft" 
regions. The hard region, in which all 3 final state particles in ([!]) are resolved, is defined such 
that S13 > s m in or S23 > s min . (In an appendix we discuss this definition when more than one 
color structure is present). The complementary region is soft. Let us review the approximations 
involved in PSS, following |nj. The 3 parton contribution to the fully differential decay can be 
written schematically as 



dT 3 = \MtfxdP$ 3 

= (\M 3 \ 2 x(l-9 s ) + \M 3 \ 2 x9 s ))xdPS 3 

= \M 3 \ 2 x (l-e 8 )dPS 3 + 6 8 x (7i(0 s )+T 2 (0 s )+T 3 (0 s )) , (4) 

where |7W 3 1 2 is the exact matrix element squared, and dPS 3 denotes the exact 3 particle phase 
space measure. Note we do not consider the effect of jet-algorithms here (they are implicit in 
the definition of the phase space). The slicing of phase space is indicated by the symbol 9 S , 
which is in the hard phase space region and 1 in the soft region. T\ is given by 

1\(0 S ) = 5|A^ 2 | 2 xrfPS soft dPS 2 

= R(9 S )\M 2 \ 2 x dPS 2 , (5) 

and represents the integral of the approximate matrix element I.M3I 2 — > S^-A^l 2 over the 
approximate phase space <iPS3 — > dPS so ft C/PS2. The resolution factor R(9 S ) is independent of 
the hard scattering and can be calculated analytically for a wide range of multiparton processes 
T 2 is given by 



T 2 (9 S ) = (\M 3 \ 2 -S\M 2 \ 2 ) xdPS 3 , (6) 
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and represents the integral over the exact 3-particle phase space phase space of the difference 
between the true matrix element and the approximate matrix element. T 3 is given by 

T 3 (6 s ) = S\M 2 \ 2 (dPS 3 -dPS 2 dPS soit ) , (7) 

and represents the difference between the integrals of the approximate matrix element over the 
true and approximate unresolved phase space. Note that T% contains the soft and collinear 
divergences needed to cancel the singularities of the virtual term, while T 2 and T3 are finite and 
vanish as the domain of support for 6 S is taken to zero. 



2.1 Matrix element 

The matrix elements for the NLO cross section for process @) are not very complicated, so we 
can be explicit. At lowest order we have 

dP 2 = ^^=Ne 2 q e 2 (8m 2 + 4s) dPS 2 

= ]--±-\M Born \ 2 dPS 2 (8) 
3 2ys 

where N is the number of colors, e q the fraction of the elementary charge e of the heavy quark, 
m its mass, s = q 2 and 

Note that in NLO approximation (A^) 2 in Eq. (|5|) is lA^Boml 2 - At 0(a s ) there are virtual 
and real emission contributions. The PSS method separates the latter into hard and soft 
contributions. The (spin and color-summed) matrix element for the real emission process (§) 
is 

\M 3 \ 2 = 16 e 2 q e 2 g 2 NC F I R (10) 
with g s the strong coupling, Cp = (N 2 — 1)/2N and 



m 2 s 23 m 2 s 12 4m 4 Am 2 s 12 s 2 12 s 23 s 12 

lR — % 2 2 1 1 1" n 1 

S13 Sf 3 S( 3 S 13 S 23 S 13 S 23 2s 13 S13 

m 2 m 2 Si 3 s± 3 m 2 s\ 2 4m 4 s\ 2 m 2 



■ +7T^ 2 — + I 11 ) 

s 13 s 23 Zs 23 s 23 s 23 s 23 s 23 

In the Ti term (j^) the eikonal approximation of the exact matrix element is used. The in- 
tegral over c?PS S oft is then performed analytically and added to the virtual corrections. The 
approximated matrix element in the soft region (|5|) is 

S\M 2 \ 2 = 16 e 2 e 2 g 2 NC F I s , (12) 
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where 



m 2 si2 4m 2 si2 m 2 s\2 



Am 



6 13 



4m 

■^23 



+ 



<? 2 
6 12 



(13) 

Sl3«23 «23 s 13 s 23 s 13«23 

Note that the difference of fllT]) and Cp~3f) which enters the T 2 term (||), is finite in the limit 
S13, S23 — > 0. 

The result of integrating (|12"D over dPS so f t is given in |L0| , and when added to the virtual 
contributions, gives the following finite expression for the 2 particle 0(a s ) differential cross 
section for process (]T|) 



with 



\M\ 



soft 



and 



\M\ 



2 

virt 



1 1 



3 2^ 



l-MiL + \M\l ilt ) dPS 2 



(14) 
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(15) 



(16) 



- 4s - 16m z - — (32Li 2 (x) + 64C 2 ) 
sp 



m 



+ - (8Li 2 (x) + I6C2) + -r In 2 (x) 8 2s 



-/?lnx(6s + 8 



m + 



lnx 



-32 — In (1 - x) + 8s In (1 - x) + 4s - 8m 2 

s 



Here C e = (Ann 2 /m 2 ) e /T(l - e), = yjl - 4m 2 /s and x = (1 - (3)/(l + 0). We have written 
the divergent contributions explicitly, even though they cancel between the soft and virtual 
contribution, so that the method independent (virtual) and method dependent (soft) terms can 
be easily read off. In particular we can obtain the results within the dipole method by replacing 
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the soft contribution with the integrated dipole terms. Note the logarithmic dependence on the 
slicing parameter s mm in the finite soft contribution. 



2.2 Phase space 

The spin-summed squared matrix elements of the previous section are functions of the final 
state momenta only via the invariants si 2 , s\ 3 , s 23 . The exact 3 particle phase space 

1 d 3 pi d 3 p 2 d?p 3 



dPS, 



-5 {4 \q -p! -p 2 -p 3 ) 



(17) 



(2tt) 5 2E x 2E 2 2E 3 

may be parametrized in terms of these invariants (after integrating over all remaining variables) 

1 1 



dPS, 



4s 32tt 3 

The integration limits of S23 at fixed S13 are 
1 



ds l2 ds 13 ds 23 S(s - S12 - si 3 - s 23 - 2m ) . 



23 



2 (si 3 + m? 



-S13 («i3 -8 + 2m 2 ) ± s 13 y/sl 3 - 2s 13 s - Asm 2 + s< 



(19) 



The limits of Si 3 at fixed s 2 3 are found by exchanging the indices 13 and 23. Setting s 23 = s 23 
we find the maxima of these two invariants 



max max 

S10 — Sn 



'13 - ^23 - s ~ 2m^s . 
In the soft (eikonal) approximation, the limits for s 23 simplify to 



(20) 



±,eik 
3 23 



— 2" ( — Sl 3 (2^ 2 — sj ± S13V s 2 — 4sm 2 ) 



2m 2 

( s — 2m 2 , s 

s 13 ~ o ± TT^P 



2m? 



2m? 



(21) 



The phase space boundaries for the exact and approximate cases are given by the Dahlitz plot 
in Fig. ffl. 
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Figure 1: Dahlitz plot for 5^3 as a function of S13 for exact (eq.(|i~9|), solid) and eikonal (eq.(p|), 
dashed) phase space boundaries at m = 5 GeV and s = 400 GeV 2 . 



5 



2.3 Results for PSS 



We now show some results for the fully inclusive cross section, as well as some differential 
distributions for process ([I]). We study what effect including the Tj contributions has on the 
■Smm dependence of the results, and shall see that including all Tj removes all s min dependence. 
We use as default values s = 400 GeV 2 , and m = 5 GeV. Figure [^a shows that, for the inclusive 
cross section, not including all Tj leads to s min dependence (in fact the T 2 worsens the s min 
dependence slightly here), but including T 2 and T3 relaxes all constraints on this parameter. 
This, however, comes at the expense of potentially lower numerical accuracy, particularly for 
the differential distributions to be considered below. The inclusion of the T3 term in particular 
requires a larger number of points in the Monte Carlo integration than using T\ alone, to achieve 
a given accuracy. In practice, therefore, it is common to use only the T\ in a PSS calculation, 
with an s min value small enough for the combined T 2 + T 3 contribution to be negligible. One 
must however be careful not to choose s m i n so small that numerical inaccuracies result from the 
large opposite sign soft+virtual and real emission contributions, as illustrated by Figure |]b. 
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Figure 2: (a) The s mm /m 2 dependence of the one-loop corrections to T(s,m 2 ), when including the 
T\ (dotted), T\ + T 2 (dashed), and T\ + T 2 + T3 (solid) contributions, (b) The s mm /m 2 dependence of 
the one-loop corrections to T(s, m 2 ) for the soft+virtual (spaced dotted) and the real emission (spaced 
dashed) final state contributions as well as their sum (solid) in the T\ + T 2 + T3 approximation. 



Turning to distributions, we show in Fig. |3| the single heavy quark transverse momentum 
and rapidity distributions at a small value of s m i n = O.OOlGeV 2 , computed with T\ only. We 
see the usual Jacobian peak near the kinematic maximum of the pt spectrum. In Fig. ^ we 
plot the s m i n dependence of the one-loop contributions to dT/px at two fixed values of px, one 
halfway and the other close to the kinematic maximum. The dip in the curves is an artifact 
which arises because at that s min and for the pt given, it is no longer kinematically possible 
for the full phase space in Fig. [l] to contribute. Note that the dip disappears for the exact 
T1 + T2+ T 3 case. Similar results are shown for the heavy quark rapidity distributions in Fig. |5] 
(where we show only the positive-rapidity part of the distribution). These figures show that 
the freedom to choose s min when including all Tj persists for distributions. 
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Figure 3: Differential decay widths at Born (dotted-dashed) and NLO (solid) levels, with parameters 
s = 400 GeV 2 , m = 5 GeV, and s m i n = 0.001 GeV 2 for differential variables (a) transverse momentum 
dr/dpr, (b) rapidity dT/dy [GeV]. 
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Figure 4: a) The s m i n -dependence of the one- loop contributions to dT/px for pt = 5 GeV. We plot 
the results including the T\ (dotted), T\ + T2 (dashed) and T\ + T2 + T3 (individual points with error 
bars) terms, b) The s m i n -dependence of the one-loop contributions to dT/pT for px = 8 GeV. Labels 
as in (a). 
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Figure 5: a) The s m i n -dependence of of the one-loop contributions to dT/dy at y = 0.3. Labels as in 
Fig. ||. b) The s m i n -dependence of the one-loop contributions to dT/dy at y = 0.6. Labels as in Fig. 



3 Dipole Subtraction 



In the dipole formalism one subtracts a suitable term from the real emission part and adds 
it again to the virtual correction after having performed one phase-space integration. The 
subtraction term consists of a sum of dipoles, each of which can be viewed as an emitter- 
spectator-antenna radiating a third particle. In the case at hand, there are only two dipoles. 
In one of these the heavy quark constitutes the emitter with the antiquark being the spectator. 
The second dipole has the roles of the quark and antiquark exchanged. The matrix element to 
be subtracted from the real emission part reads: 



\M> 



Born 



1 



X 



Ul3 



2(1 - r r) - (1 - r ) 



r r^(l - r)(l - r r) 
1 - r 



(22) 



Here 



r = /3 2 



$13 + s 23 

s — Am 2 ' 



1 — UqU 
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2(1 - r r) - (1 - r ) - 



1 ~ r 
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1 s 



13 



^23 



1 — ror 



U Si3 + S 2 3 



(23) 



This contribution is then integrated over the dipole phase space and added to the virtual 
corrections. The integrated version reads: 



dPS 




1-Ii±^lni± 
2 1- 



—2 lnro — In 



S 



' f r ° * Li 2 (V^) - Li 2 (-V^) + 2 Li 2 ( l±^) -2UJ 1 '^ 



_21nr ln f 1 + v / ^ V ln21n v 7 ^ + 1^2 

u-^y i + v^ 2 



+ ln(l - Vfb) ^ [ — J + - ln 2 (l + 0^) - ^ ln 2 (l - v^) 



+0(e) (24) 

The poles in e cancel against those of the virtual corrections. 

We do not show separate results for the dipole method, which is exact and independent of 
any theoretical cut-off parameter. Numerical results for the dipole method in comparison to 
the PSS method can be found in the next section. 



4 Comparisons of PSS and dipole subtraction 

In this section we perform some numerical comparisons between the two methods for the process 



at hand. We use as phase space measure the expression (|18|). The integrations over its variables 
are performed using the well-known Monte Carlo iterative integration routine VEGAS [20] . We 
note that we found similar results when we used (|TTD, generating the 4- vectors via a cascade 
algorithm. This required using more random number points in order to achieve the same 
accuracy. 

The PSS method is relatively easy to implement, with little analytical calculation, at the 
expense of requiring cancellations between large numbers (for small s m j n ) or having multiple 
negative contributions (for large s min when including T 1; T 2 and T 3 ). Since the dipole method 
requires more analytical preparation work to be implemented, we expect it to show better 
numerical integration in the Monte Carlo program. We will see that this expectation is borne 
out by our results. 

Our first comparison addresses the relative accuracy achieved in the computation of the 
inclusive cross section as a function of the number of points, for 20 iterations, of which we 
use the first five to set the VEGAS grid pOj, leaving a sample of N = 15 results. For each 
method, we perform separate runs for the 0(a s ) 2-particle and 0(a s ) 3-particle contributions, 
and combine them for each iteration, leading to 15 results r^. The mean result r and its error 
5r are then computed as 



1 N 



r 



, \ <y> . Of \ • 

N fr[ *' \N N 



i Ef =1 (n - rf 



(25) 



The results for this comparison are given in Table [1]. We note that the PSS method suffers 
further penalties in accuracy and efficiency if the value of s m i n is chosen so large that the T 2 
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s = 400 GeV 2 


points 


DIP 


PSS 


1000 


0.04% 


1% 


10000 


0.009% 


0.3% 


100000 


0.003% 


0.1% 



Table 1: Accuracy 5r/r of the inclusive cross section attained for a given number of points per 
iteration in the two methods. The same phase space and random number generators are employed. 
The PSS results use the T\ contribution only, with s m j n = O.OOlGeV 2 . 



iteration 


DIP 


PSS 


1 


0.1% 


100% 


2 


0.09% 


70% 


3 


0.06% 


10% 



Table 2: Comparison of the two methods as to the approximate relative deviations of their first three 
(grid-setting) iterations from the final mean (computed starting from the fifth iteration), for the case 
of the inclusive cross section. The same phase space and random number generators are employed, at 
s = 400 GeV 2 . The PSS results use the T x contribution only, with s min = O.OOlGeV 2 . 



and T 3 become necessary; in particular the T 3 contribution requires generating the soft phase 
space measure, and involves the difference of two phase space measures which are very similar 
in magnitude for small values of the soft invariants, cf. Eq. (^). 

Our second comparison addresses the efficiency in the computation of the inclusive cross 
section as a function of the number of iterations, for 10 4 random number points. We see that 
the dipole method reaches a given accuracy with less iterations. 

Next we compare the efficiency of these methods to compute transverse momentum and 
rapidity distributions. As before, the first 5 of 20 iterations are used solely for grid-setting, 
with 10 4 points per iteration. The values and their errors for each bin are computed according 
to (|25|) . We see that the dipole methods produces somewhat smaller errors, with slightly less 
bin to bin variations. When we increased the number of points, we saw that both methods 
perform not too differently. We also noticed that this loss of superiority is progressive with the 
number of bins. This suggests that parts of the positive and negative contributions end up in 
different bins. To test this idea, we performed a simple smearing where each event with weight 
w that would normally end up in bin % is distributed in bins % — 1, i, % + 1, each with weight w/3. 
We found that this smearing indeed reduced the errors somewhat, but in about equal measure 
for both methods. 

Finally, we compared the accuracies of the methods in the large s limit for the inclusive 
cross section, where the heavy quarks become effectively massless. For the PSS method we 
investigated how to choose s m i n in order to minimize the calculation error. We found that 
s min is best chosen not too small (which would lead to large numerical cancellations), and as a 
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Figure 6: Differential decay widths at NLO level, with parameters s = 400 GeV 2 , m = 5 GeV, for the 
phase space slicing method (at s min = O.OOlGeV 2 ) for differential variables (a) transverse momentum 
dT/dpT (b) rapidity dT/dy [GeV] 
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7: Differential decay widths at NLO level, with parameters s = 400 GeV 2 , m = 5 GeV, 
dipole method for differential variables (a) transverse momentum dT/dpx (b) rapidity dT/dy 
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fraction of s between 0.01 and 0.1. For large s m j n this may require the inclusion of the T 2 and 
T 3 terms, which contribute about 10% to the cross section for s m j n = 0.1 s at s — 250000 GeV 2 . 
For s = 250000 GeV 2 and s m i n = 0.01 s their contribution is only 2%, with a slightly larger 
total error for the same number of points. We found similar results keeping s fixed and letting 
m become smaller. For the dipole method we found that it has consistently better accuracy 
than the PSS method in these limits. In general in the high-energy limit, both these methods 
lead to cancellations between contributions with + ln(s/m 2 ) and — ln(s/m 2 ) terms, which is not 
advantageous numerically. Therefore a method which avoids such logarithms could be desirable 



5 Conclusions 

In this paper we have compared the accuracy and efficiency of two general-purpose methods 
to compute NLO heavy quark production cross sections, for a very simple case. We found the 



dipole method JT2L 13], while involving additional analytical work, to be superior in efficiency 



and accuracy. A similar conclusion was reached by Dittmaier |L4j who compared his method 
with a slicing calculation for a number of electroweak cross sections. 

The phase space slicing method || 0, 110], which is easy to use and minimizes analytical 



work, can be extended ||T9[ to become fully independent of the slicing parameter, which we 
demonstrated in this paper for the reaction at hand. Although our case-study involves only the 
simplest of heavy quark production processes, it is, we believe a useful first step toward gaining 
numerical experience with general methods for constructing NLO Monte Carlo programs for 
heavy quark production. Moreover, such experience gained at NLO is likely to be very valuable 
when these methods are generalized for NNLO cross sections. 
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A A subtlety in Phase Space Slicing 

In this appendix we mention a subtlety in implementing PSS, which is well known to experts in 
the field, but not readily found in the literature. It does not come into play for our simple case 
study, but does for amplitudes in which more than one color structure is present. The correct 
implementation of phase space slicing requires that the real emission amplitude be decomposed 
into pieces with a unique singularity structure and that the slicing procedure be defined for 
each piece separately. 

We discuss the point for a simplified example, where we take al particles to be massless. 
We assume that the real emission amplitude is given by 

\M,\ 2 = + . (A.l) 
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This simple example has the singularity structure for the leading-color part of the 7* — > qqgg 
amplitude. The correct way to implement phase space slicing treats each singularity structure 
separately and the resolved contribution from the real emission amplitude reads therefore 



/ d PS 4 9 (SX3 — S m j n ) 9 (s 34 — S mm ) 9 (s 2 4 — Smin) 

J S13S34S24 

+ / c/PS 4 9 (si4 - s min ) 9 (s 34 - s min ) 9 (s 23 - s min ) . (A.2) 

J S14S34S23 

The incorrect method cuts out all possible singularities from the amplitude and uses the ex- 



pression 



/ 



c/PS 4 ( 1 ) 9 (s 13 - s min ) 6 (s 34 - s min ) 9 (s 24 - s r 



S13S34S24 Si 4 S 34 S 2 3 



• 9 (su - s min ) 9 (s 23 - s min ) (A. 3) 



for the resolved contribution. The difference between the correct implementation and the 
incorrect one consists of terms of the form 

f d PS 4 9 - s min ) 9 (s 34 - s min ) 9 O24 - Smin) 9 (s min - s u ) 9 (s 23 - s u ) (A.4) 

J S13S34S24 

together with three similar terms, obtained by exchanging 3 4 in the matrix element and 
the slicing procedure. Contrary to naive expectations, these terms do not vanish in the limit 
Smin — * 0, but give a constant contribution. 
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